#include "Analysis.hh"
#include "TMath.h"
#include "RawData.hh"
#include "HeaderParameters.hh"
#include <string>
#include <sstream>

Analysis::Analysis(string rootfilename, HeaderParameters* header_arg, RawData* my_data_arg, float voltthreshold, float varthreshold, int channels){

  rootfile = new TFile(rootfilename.c_str(), "RECREATE");
  tree = new TTree("tree", "tree");

  char chnum[100];

  char chargename[100];
  //~ char noisename[100];
  //~ char variancename[100];
  char vmaxname[100];
 
  char chargevarname[100];
  char vmaxvarname[100];
  //~ char noisevarname[100];
  //~ char variancevarname[100];

  int ch;
  for(ch=0;ch<channels;ch++)
  {
	  strcpy(chargename, "charge");
	  strcpy(vmaxname, "vmax");
	  //~ strcpy(noisename, "noise");
	  //~ strcpy(variancename, "variance");

	  sprintf(chnum, "%d", ch);

	  strcat(chargename, chnum);
	  strcat(vmaxname, chnum);
	  //~ strcat(noisename, chnum);
	  //~ strcat(variancename, chnum);

	  strcpy(chargevarname, chargename);
	  strcpy(vmaxvarname, vmaxname);
	  //~ strcpy(noisevarname, noisename);
	  //~ strcpy(variancevarname, variancename);

	  strcat(chargevarname,"/F");
	  strcat(vmaxvarname,"/F");
	  //~ strcat(noisevarname,"/F");
	  //~ strcat(variancevarname,"/F");
	  
	  tree->Branch(chargename, &charge[ch], chargevarname);
	  tree->Branch(vmaxname, &trace_vmax[ch], vmaxvarname);
	  //~ tree->Branch(noisename, &noise[ch], noisevarname);
	  //~ tree->Branch(variancename, &variance[ch], variancevarname);
  }

	
	
  //-------general use---------//
//  time_per_div = header->time_per_div;
//  num_traces = header->num_traces;
//  num_traces_cut = num_traces;
//  num_bins = int(header->GetParameter(3));
	my_channels=channels;
	trace_electrons = new TH1F*[2];
	voltage_noise = new TH1F*[2];
	voltage_maxnoise = new TH1F*[2];
	header = header_arg;
  voltage_threshold = voltthreshold;
  variance_threshold = varthreshold;
  my_data = my_data_arg;
  TVectorF vthresh(1);
  vthresh[0] = voltage_threshold;
  vthresh.Write("vthresh");
  //---------------------------//
  //--------noise histo---------//
  first_sample_noise=true;
  total_v_noise = 0;
  num_over=0;
  //---------------------------//
  //------variance histo-------//
  first_trace_variance = true;
  total_sigma = 0;
  //---------------------------//
  //--------count peaks--------//
  num_peaks=0; 

  for(ch=0;ch<my_channels;ch++)
  {
	num_non_empty[ch]=0;
  }
  //---------------------------//
  //-------charge_histo--------//
  first_trace_charge=true;
  num_neg=0;
  total_num_electrons[0] = 0;
  total_num_electrons[1] = 0;
  all_electrons_cut = new bool[num_traces];
  num_electrons_threshold = -7.981e8;
  num_electrons_outside = 0;
  //---------------------------//


}

Analysis::~Analysis(){
  //~ tree->Write();
  //~ rootfile->Close(); 
  rootfile->Write();
  delete [] all_electrons_cut;
  delete [] trace_electrons;

}



//-------charge_histo-------//

//void Analysis::charge_minmax()
//{
	//my_waveforms = my_data->GetCurrentWaveforms();
	//for(int channel=0;channel<my_channels;channel++)
	//{
		//float time_per_div = header->GetParameter(channel,5);
		//charge = float(my_waveforms[channel]->Integral())/50.0*time_per_div;
	  	//charge*=1E9; //Convert to Nanocoulombs
	
		////num_electrons= charge/1.6e-19; //1.6e-19 is the charge of an electron
		//num_electrons= charge; //1.6e-19 is the charge of an electron
	
		//if (first_trace_charge == true)
		//{
	    //max_num_electrons[channel] = num_electrons;
	    //min_num_electrons[channel] = num_electrons;
		//first_trace_charge = false;
		//}
	
	    //if (num_electrons > max_num_electrons[channel])
	    //{
	    //max_num_electrons[channel] = num_electrons;
	    //}
	
	    //if (num_electrons <  min_num_electrons[channel])
	    //{
	    //min_num_electrons[channel] = num_electrons;
	    //}
	
		//total_num_electrons[channel]+=num_electrons; 
	//}

	//max_total_electrons = 0;
	//for (int channel=0; channel<my_channels; channel++)
	    //max_total_electrons += max_num_electrons[channel];
//}


//void Analysis::charge_range()
//{
	//std::stringstream chan_num;
	//std::string hist_name;
	//for(int channel=0;channel<my_channels;channel++)
	//{
	////	float maxrange = pow(10,ceil(log10(fabs(max_num_electrons)))) * (max_num_electrons/fabs(max_num_electrons));
	////	float minrange = pow(10,floor(log10(fabs(min_num_electrons))))* (min_num_electrons/fabs(min_num_electrons));

		//int num_traces = int(header->GetParameter(channel,7));
	
		//mean_num_electrons = total_num_electrons[channel]/num_traces;
	
		//float maxrange = mean_num_electrons + 2*fabs(max_num_electrons[channel]-mean_num_electrons);
		//float minrange = mean_num_electrons - 2*fabs(min_num_electrons[channel]-mean_num_electrons);

		//chan_num << channel;
		//hist_name = "trace_electrons"+chan_num.str();
		//chan_num.str("");
	    ////trace_electrons[channel] = new TH1F(hist_name.c_str(), "Charge per event", 100, minrange, maxrange);
	    //trace_electrons[channel] = new TH1F(hist_name.c_str(), "Charge per event", 100, minrange, maxrange);
	//}
	//total_trace_electrons = new TH1F("summed_charge", "Quantum Dot Nanocrystals with Po210", 100, 0, max_total_electrons);
	//total_trace_electrons->GetXaxis()->SetTitle("Nanocoulombs");
	//total_trace_electrons->GetYaxis()->SetTitle("Number of Events");
	//total_trace_electrons->GetXaxis()->CenterTitle();
	//total_trace_electrons->GetYaxis()->CenterTitle();
//}

//void Analysis::charge_histo(bool writewaveform){
	  //my_waveforms = my_data->GetCurrentWaveforms();
	//float totalcharge = 0;
	//for(int channel=0;channel<my_channels;channel++)
	//{
		
	  //float time_per_div = header->GetParameter(channel,5);
	  //charge = float(my_waveforms[channel]->Integral())/50.0*time_per_div; //Coulombs, I=dq/dt V=(dq/dt)*R integral(V) dt / R = q, R from the oscilloscope = 50 ohms, division by timescale parameter to convert current from coulombs/timediv to coulombs/sec.
	  //charge*=1E9; //Convert to Nanocoulombs
	 //totalcharge += charge;
	 //// cout  << "NENT: " << my_waveforms[channel]->GetEntries() << endl;
	 //// cout  << "INT: " << my_waveforms[channel]->Integral() << endl;

	 //// cout << "CHARGE" << charge << endl;
	
	  ////num_electrons= charge/1.6e-19; //1.6e-19 is the charge of an electron
	  //num_electrons= charge; //1.6e-19 is the charge of an electron
	
	  //if (num_electrons < num_electrons_threshold)
	  //{
	  	//num_electrons_outside++;
	  //}
	
	  //trace_electrons[channel]->Fill(num_electrons); // Fill the my_waveforms electron histogram
	  
	  //if (writewaveform)
	      //my_waveforms[channel]->Write(); // Write the waveform to the root file

		  ////my_waveforms->Smooth();
	      ////my_waveforms->Write(); // Write the waveform to the root file
	
	//// ----Look at negative my_waveforms---
	//
	      //if (num_electrons < -5e27)
	      //{
	////    writemy_waveforms = true;
	////    analysis->AnalyzeWaveform(my_data, time_per_div, &num_electrons, writemy_waveform);
	
	      //num_neg++;
	////    my_waveforms->SaveAs("
	      //}
	//*/
	////-------------------------
	//}

	//total_trace_electrons->Fill(totalcharge);


//}

void Analysis::charge_histo(bool writewaveform)
{
	my_waveforms = my_data->GetCurrentWaveforms();

	for(int channel=0;channel<my_channels;channel++)	
	{		
		float time_per_div = header->GetParameter(channel,5);
		
		charge[channel] = float(my_waveforms[channel]->Integral())/50.0*time_per_div; //Coulombs, I=dq/dt V=(dq/dt)*R integral(V) dt / R = q, R from the oscilloscope = 50 ohms, division by timescale parameter to convert current from coulombs/timediv to coulombs/sec.
		
		charge[channel]*=1E9; //Convert to Nanocoulombs
		
		if (writewaveform)
		my_waveforms[channel]->Write(); // Write the waveform to the root file	
	}
}


//~ void Analysis::write_charge_histo()
//~ { 
	//~ for(int channel=0;channel<my_channels;channel++)
	//~ {
		//~ trace_electrons[channel]->Write();
		//~ cout<< "---------------charge_histo-----------------"<<endl;
	    //~ cout<< max_num_electrons[channel] << " " << min_num_electrons[channel] << endl;
	    //~ cout<< "max num electrons" << max_num_electrons[channel] << endl;
	    //~ cout<< "min num electrons" << min_num_electrons[channel] << endl;
		//~ cout<< "Num Neg: " << num_neg << endl;
		//~ cout<< "Noise threshold (3 sigma below mean): " << num_electrons_threshold << endl;
		//~ cout<< "Number of traces below threshold: " << num_electrons_outside << endl;
		//~ cout<< "Percentage of traces below threshold: " << 100*(float(num_electrons_outside)/float(num_traces)) << " %" << endl;
		//~ cout<< "------------------------------------------"<<endl;
		//~ cout<< endl;
	//~ }
		//~ total_trace_electrons->Write();
//~ }

//--------------------------//
////*
//-------noise histo--------//

//~ void Analysis::noise_minmax()
//~ {
	//~ for(int channel=0;channel<my_channels;channel++)
	//~ {
		//~ my_waveforms = my_data->GetCurrentWaveforms();
		//~ int num_bins = int(header->GetParameter(channel,3));
	//~ 
		//~ if (first_sample_noise == true)
		//~ {
			//~ voltage = my_waveforms[channel]->GetBinContent(1); 
			//~ min_v_noise = voltage;
			//~ max_v_noise = voltage;
			//~ first_sample_noise = false;
		//~ }
	 //~ 
		//~ for(int count=1; count<=num_bins; count++)
		//~ {
			//~ voltage = my_waveforms[channel]->GetBinContent(count); 
	//~ 
	        //~ if (fabs(voltage) > voltage_threshold)
	        //~ {
	        	//~ num_over++;
	        //~ }
	//~ 
	        //~ if  (voltage < min_v_noise)
	        //~ {
	            //~ min_v_noise = voltage;
	        //~ }
	//~ 
	        //~ if  (voltage > max_v_noise)
	        //~ {
	            //~ max_v_noise = voltage;
	        //~ }
			//~ total_v_noise += voltage;
		//~ }
	//~ }
//~ }
//~ 
//~ void Analysis::noise_range()
//~ {
	//~ for(int channel=0;channel<my_channels;channel++)
	//~ {
		//~ int num_traces = int(header->GetParameter(channel,7));
		//~ int num_bins = int(header->GetParameter(channel,3));
		//~ int num_traces_cut = num_traces;
		//~ mean_v_noise = total_v_noise/(num_traces_cut*num_bins);
	//~ 
		//~ float maxrange = mean_v_noise + 2*fabs(max_v_noise-mean_v_noise);
		//~ float minrange = mean_v_noise - 2*fabs(min_v_noise-mean_v_noise);
//~ //		float maxrange = max_v_noise;
//~ //		float minrange = min_v_noise;	
//~ 
		//~ float absmax;
//~ 
		//~ if (fabs(max_v_noise)> fabs(min_v_noise))
		//~ {
			//~ absmax = fabs(max_v_noise);
		//~ }
		//~ else
		//~ {
			//~ absmax = fabs(min_v_noise);
		//~ }
//~ 
		//~ voltage_noise[channel] = new TH1F("voltage_noise", "Voltage noise", 100, minrange, maxrange);
		//~ voltage_maxnoise[channel] = new TH1F("voltage_maxnoise", "Max |voltage| per trace",100, 0, absmax);
	//~ }
//~ }

//~ void Analysis::noise_histo()
//~ {
	//~ float max_volts;
//~ 
	//~ for(int channel=0;channel<my_channels;channel++)
	//~ {
		//~ my_waveforms = my_data->GetCurrentWaveforms();
		//~ int num_bins = int(header->GetParameter(channel,3));
		//~ if (first_sample_noise == true)
		//~ {
			//~ voltage = my_waveforms[channel]->GetBinContent(1); 
	//~ //		min_v_noise = voltage;
	//~ //		max_v_noise = voltage;
			//~ max_volts=fabs(voltage);
			//~ first_sample_noise = false;
		//~ }
	    //~ for(int count=1; count<=num_bins; count++)
	    //~ {
	//~ 
			//~ voltage = my_waveforms[channel]->GetBinContent(count); 
	//~ 
	        //~ voltage_noise[channel]->Fill(voltage);
	//~ /*
	        //~ if (voltage > voltage_threshold)
	        //~ {
	        	//~ num_over++;
	        //~ }
	//~ 
	//~ */
	        //~ if  (fabs(voltage) > max_volts)
	        //~ {
	            //~ max_volts = fabs(voltage);
	        //~ }
	    //~ }
//~ 
		//~ voltage_maxnoise[channel]->Fill(max_volts);
//~ 
		//~ }
	//~ 
//~ }

//~ void Analysis::write_noise_histo()
//~ {
	//~ for(int channel=0;channel<my_channels;channel++)
	//~ {
		//~ voltage_noise[channel]->Write();
		//~ voltage_maxnoise[channel]->Write();
	//~ 
		//~ cout<< "--------------voltage noise---------------"<<endl;
		//~ cout<< "threshold: " << voltage_threshold << " Volts"<< endl;
		//~ cout<< "Number of samples over threshold: " << num_over << endl;
		//~ cout<< "Largest voltage: " << max_v_noise << endl;
		//~ cout<< "Smallest voltage: " << min_v_noise << endl;
		//~ cout<< "------------------------------------------"<<endl;
		//~ cout<< endl;
	//~ }
//~ }


//--------------------------//
////---------variance_cut---------//
//
//bool Analysis::variance_cut(bool first_loop)
//{
//	my_waveforms = my_data->GetCurrentWaveforms();
//	for(int channel=0;channel<my_channels;channel++)
//	{
//	
//		v_total = 0;
//		v_mean = 0;
//		v_diff_total = 0;
//		sigma_squared = 0;
//	
//	    for(int count=0; count<(num_bins/10); count++)
//	    {
//			voltage = my_waveforms[channel]->GetBinContent(count); 
//	
//			v_total += voltage;
//	    }
//	
//		v_mean = v_total/(num_bins/10);
//	
//	    for(int count=0; count<(num_bins/10); count++)
//	    {
//			voltage = my_waveforms[channel]->GetBinContent(count); 
//	
//			v_diff_total += pow(voltage - v_mean, 2);
//	    }
//	
//		sigma_squared = v_diff_total/(num_bins/10);
//		sigma = sqrt(sigma_squared);
//	
//		if(sigma > variance_threshold)
//		{
//			if(first_loop==false)
//			{
//			num_traces_cut --;
//			}
//			return false;
//		}
//		
//		return true;
//	}
//}
//
//void Analysis::write_variance_cut()
//{
//	for(int channel=0;channel<my_channels;channel++)
//	{
//		cout << "-----------variance_cut-----------"<<endl;
//		cout << "Number of traces: " << num_traces << endl;
//		cout << "Number of traces cut: " << num_traces - num_traces_cut << endl;
//		cout << "Number of traces left: " << num_traces_cut << endl;
//		cout << "Percentage cut: " << 100*float(num_traces - num_traces_cut)/float(num_traces) << " %" << endl;
//		cout << "----------------------------------"<<endl;
//		cout << endl;
//	}
//}
//
////--------------------------//
//
////---------variance_histo---------//
//
//void Analysis::variance_minmax()
//{
//	for(int channel=0;channel<my_channels;channel++)
//	{
//		my_waveforms = my_data->GetCurrentWaveforms();
//	
//		v_total = 0;
//		v_mean = 0;
//		v_diff_total = 0;
//		sigma_squared = 0;
//	
//	    for(int count=1; count<=(num_bins/10); count++)
//	    {
//			voltage = my_waveforms[channel]->GetBinContent(count); 
//	
//			v_total += voltage;
//	    }
//	
//		v_mean = v_total/(num_bins/10);
//	
//	    for(int count=1; count<=(num_bins/10); count++)
//	    {
//			voltage = my_waveforms[channel]->GetBinContent(count); 
//	
//			v_diff_total += pow(voltage - v_mean, 2);
//	    }
//	
//		sigma_squared = v_diff_total/(num_bins/10);
//		sigma = sqrt(sigma_squared);
//	
//		if (first_trace_variance == true)
//		{
//			max_sigma = sigma;
//			min_sigma = sigma;
//			first_trace_variance = false;
//		}	
//	
//		if (sigma > max_sigma)
//		{
//			max_sigma = sigma;
//		}
//	
//		if (sigma < min_sigma)
//		{
//			min_sigma = sigma;
//		}
//	
//		total_sigma += sigma;
//	}
//}
//
//void Analysis::variance_range()
//{
//	for(int channel=0;channel<my_channels;channel++)
//	{
//		mean_sigma = total_sigma/num_traces;
//	
//		float maxrange = mean_sigma + 2*fabs(max_sigma-mean_sigma);
//		float minrange = mean_sigma - 2*fabs(min_sigma-mean_sigma);
//	
//		variance = new TH1F("variance", "Variance", 100, minrange, maxrange);
//	}
//}
//
//void Analysis::variance_histo()
//{
//	for(int channel=0;channel<my_channels;channel++)
//	{
//	
//		my_waveforms = my_data->GetCurrentWaveforms();
//	
//		v_total = 0;
//		v_mean = 0;
//		v_diff_total = 0;
//		sigma_squared = 0;
//	
//	    for(int count=1; count<=(num_bins/10); count++)
//	    {
//			voltage = my_waveforms[channel]->GetBinContent(count); 
//	
//			v_total += voltage;
//	    }
//	
//		v_mean = v_total/(num_bins/10);
//	
//	    for(int count=1; count<=(num_bins/10); count++)
//	    {
//			voltage = my_waveforms[channel]->GetBinContent(count); 
//	
//			v_diff_total += pow(voltage - v_mean, 2);
//	    }
//	
//		sigma_squared = v_diff_total/(num_bins/10);
//		sigma = sqrt(sigma_squared);
//	
//	//	cout << "sigma squared: " << sigma_squared << endl;
//	
//		variance ->Fill(sigma);
//	}
//}
//
//void Analysis::write_variance_histo()
//{
//	for(int channel=0;channel<my_channels;channel++)
//	{
//		variance->Write();
//	
//		cout<< "-----------------variance-----------------"<<endl;
//	//	cout<< "threshold: " << variance_threshold << " Sigma"<< endl;
//	//	cout<< "Number of samples over threshold: " << num_over << endl;
//	//	cout<< "Largest sigma: " << max_v << endl;
//	//	cout<< "Smallest sigma: " << min_v << endl;
//		cout<< "------------------------------------------"<<endl;
//		cout<< endl;
//	}
//}



//---------count peaks------//

bool Analysis::count_peaks(bool first_loop_peakcount)
{
	my_waveforms = my_data->GetCurrentWaveforms();
	for(int channel=0;channel<my_channels;channel++)
	{
		trace_vmax[channel]=my_waveforms[channel]->GetMaximum();
		if(trace_vmax[channel]>=voltage_threshold)
		{
			num_non_empty[channel]++;
		} 
//		int num_bins = int(header->GetParameter(channel,3));
//		float voltage;
//	//	bool in_peak = false;
//	
//	
//		//my_waveforms->Smooth(5); //Smooth the waveform to reduce peak triggers by noise
//	
//	
//	
//	    for(int count=1; count<=num_bins; count++)
//	    {
//	
//			voltage = my_waveforms[channel]->GetBinContent(count); 
//		
//
//
//	
//			if (fabs(voltage) > voltage_threshold)
//			{
//				if (first_loop_peakcount == true)
//				{
//					num_non_empty++;
//				}
//				//break;
//				return false;
//			}
//	
////			if (voltage > threshold && in_peak == false)
////			{
////				in_peak = true;
////			}
////	
////			if (voltage <= threshold && in_peak == true)
////			{
////				in_peak = false;
////				num_peaks++;
////			}
//				
//	    }
//	
//		return true;
		return true;
	}

}	

void Analysis::write_count_peaks()
{

	for(int channel=0;channel<my_channels;channel++)
	{
		int num_traces = int(header->GetParameter(channel,7));
		int num_traces_cut = num_traces;
		cout << endl;
		cout <<"@@@@---Channel: "<<channel<<"---@@@@"<<endl;
		cout<< "---------------peak count-----------------"<<endl;
		cout<< "Threshold: " << voltage_threshold << " Volts"<< endl;
		//cout<< "Number of peaks: " << num_peaks << endl;
		cout<< "Number of empty traces: " << num_traces_cut - num_non_empty[channel] << endl;
		cout<< "Number of traces (after cut): " << num_traces_cut << endl;
		cout<< "Percentage empty: " << (float(num_traces_cut - num_non_empty[channel])/float(num_traces_cut))*100 << " %" << endl;
		cout<< "------------------------------------------"<<endl;
		cout<< endl;
		cout << endl;
	}
}

//--------------------------//

void Analysis::filltree(void)
{
	tree->Fill();
}

void Analysis::WriteRootfile(){
 // rootfile->Write();
 // my_waveforms->Write();
}


//void Analysis::PlantTree(){

  //tree->Branch("Qall", &Qall, "Qall/F");
  //tree->Branch("Qpeak", &Qpeak, "Qpeak/F");


//}


//void Analysis::ClearRootVariables(){
//  Qall = 0;
//  Qpeak = 0;
//}
